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SUMMARY 


A new eigensolution routine, FEER (Fast Eigensolution Extraction 
Routine) , used in conjunction with NASTRAN at Israel Aircraft Industries 
is described. The FEER program is based on an automatic matrix re- 
duction scheme whereby the lower modes of structures with many degrees 
of freedom can be accurately extracted from a tridiagonal eigenvalue 
problem whose size is of the same order of magnitude as the number of 
required modes. The process is effected without arbitrary lumping of 
masses at selected node points or selection of nodes to be retained in 
the analysis set. 

The results of computational efficiency studies are presented, 
showing major arithmetic operation counts and actual computer run times 
of FEER as compared to other methods of eigenvalue extraction, including 
those available in the NASTRAN READ module. It is concluded that the 
tridiagonal reduction method used in FEER would serve as a valuable addi- 
tion to NASTRAN for highly increased efficiency in obtaining structural 

vibration modes. 


INTRODUCTION 


One of the most burdensome computational tasks in discretized 
structural systems centers around the extr ution of mode shapes and 
frequencies when the orders of the matr^' -a are large. The difficulties 
are compounded as the number of required eigensolutions increases and 
multiple or near-multiple eigenvalues are encountered. 
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Currently, NASTRAN provides three methods for modal extraction 
(refs. 1 and 2): the Tridiagonal or Givens method, the Inverse Power 

Method with Shifts, and the Determinant method. In each method the 
problem size encountered is equal to the number of degrees-of-freedom 
in the analysis set which, given typical, present-day problem applica- 
tions, may run into the thousands. One means of reducing the size of 
the analysis set is via the Guyan reduction (ref. 3), which has been 
incorporated into NASTRAN. This technique, which is similar in concept 
to the Kaufman - Hall reduction (ref. 4), requires a "judicious" 
elimination of selected mass degrees-of-freedon. and an attempt is made 
to account for the influence of the eliminated nodes through equivalent 
energy criteria. As demonstrated by Levy (ref. 5), such an intuitive 
approach involves a great deal of guesswork and can le d to grossly 
inaccurate results, particularly in systems with relatively non-uniform 
mass distribution; . 


What is required to circumvent these difficulties is a more 
automated eigenreduction scheme which yields accurate lower modes of 
the structural system. In essence, the probli m may be post^ as follows: 


Given the nth order eigenvalue problem 

[K] {$} - w 2 [M){*} (1) 

where [K] and [M] are symmetric and non-negative definite, we wish to 
approximate the modal vectors by 

4) • (T]{«) (2) 


where [T] is a suit*Mv constructed transformation matrix of size 
nxm (m«n) and {6 1 is f.n ra- component vector of generalized co- 
ordinates. Usin; a Rhyj eigh-Ritz procedure the resulting reduced, 
mth order eigenpt'iblem is of the form 


where 


and w 
lover 
order 
rich 



» u 2 [*H6) 

(3) 

CKJ - 

[t] t ;k] m 

(4a) 

[MJ - 

[T] 7 [Ml IT] 

(4b) 

is an 

approximate modal frequency. 

If a specified number of 


modes are to be eccuretely obtained, then the indivlduel nth 
vectors comprising the trensformation matrix must be sufficiently 
in the corresponding modal vecto r .. Thus, the practical value of 
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the reduction 3cheme hinges on its ability to generate such a trans- 
formation matrix with a minimum of corputational effort. 

A number of closely related methods involving eigenreduction 
concepts have been proposed previously. In the work of Hestenes and 
Karush (ref. 6), eigensolutions were obtained via a block power method 
(Iterating with several vectors simultaneously as opposed to a single 
vector) and a reduced eigenvalue problem was employed to orthogonalize 
and Improve successive blocks of vectors between iteration steps. More 
recently, Jennings and Orr (ref. 7), Dong, Wolf and Peterson (ref. 8), 
and Bathe and Wilson (ref. 9) proposed similar techniques using the 
Inverse Power Method in conjunction with simultaneous sets of vectors 
(alternately called Simultaneous Iteration, Subspace Iteration and 
Block-Stodola methods). In each of these approaches, however, the 
functional role of the reduced eigenproblem is to improve a sub- 
space of approximate modal vectors with central emphasis being placed 
on a block-type Inverse Power method. 

In this report, a new eigenreduction routine, FEER (Fast 
Eigensolution Extraction Routine) is described, wherein a sing 1 * reduced 
eigenproblem is generated for the accurate extraction of any specified 
number of lower modes. Further, the transformation matrix is generated 
vector-by-vector in such a way that the reduced eigenproblem is 
tridiagonal in form. The FEER program is now being used in conjunction 
with NASTRAN at Israel Aircraft Industries to obtain much more 
economical eigensolutions than currently possible with the NASTRAN READ 
module. 


The tridiagonal reduction method employed in FEER was first 
suggested by Crandall (ref. 10) as a truncated version of the Lanczos 
algorithm (ref. 11). However, it was soon discovered that the original 
scheme possessed numerical instabilities (refs. 12 and 13). The 
necessary Improvements to correct these weaknesses were made by Ojalvo 
and Newman (ref. 14), who were the first to develop a successful 
tridiagonal reduction prog' a for large scale structural dynamics 
problems. The FEER computer program contains further refinements later 
Introduced by Newman and Pipano (ref. 15), including: 

1. Highly efficient numerical computation schemes, uslvg packing 
techniques which take advantage of matrix sparsity. 

2. Calculation of accurate upper and lower bounds on the 
extracted eigenvalues. 

3. Accommodation of singular mass matrices and stiffness matrix 
singularities associated with rigid body modes. 
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TRIDIAGONAL REDUCTION METHOD 


Preliminary Operations 


Employing the NASTRAN notation, the structural eigenvalue problem 
i9 of the form 




( 5 ) 


Both [K ] and [M 1 are nth order symmetric , non-negative and seni- 
aa aa 


definite untrices corresponding to the analysis set. Hence, they may 
both be singular, but all vhe eigemal es are zero or positive. 


In order to obtain a decomposable matrix, a small, positive shift 
parameter, a , is chosen such that 


0) 


2 

o - o 

o 


( 6 ) 


Then 


[K ^ ** H* J * 

aa aa a o aa a 


( 7 ) 


It can be easily show . that the shifted stiffness matrix is non- 
singular and positive-definite provided that ihe system masses generate 
kinetic energy due to any kinematically admissible rigid body motions 
of the structure. This requirement is always satisfied by the mass 
matrix in a physically veil-posed problem. 


In order to maintain the elements of the subsequent trial vectors 
on the order of unity, a positive mass-scaling paraiuSter, S, in also 
employed, such that 


lM aa J .-f-lM] 


( 8 ) 


If a Cholasky symmetric decomposition of the shifted stiffness 
matrix is performed: 




in nr 

it follows that the eigenvalue problem, squat ion (7). be converted 


(9) 
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to the form 


[B] {x} = X{x} (10) 

where 

[B] = [L]‘ 1 [M] [L _1 ] T (lla) 

{x} = [L] T {* } (lib) 

a 

and 

X = ~~ (lie) 

U) +0 

The above triangular matrix inverses are treated as purely 
operational symbols, since in actual numerical operations forward and 
backward passes on vectors are employed. 

Generation of the Reduced Eigenproblem 


A reduction of the nth order eigenvalue problem, equation (10), is 
effected through the transformation 

{ x } = [ V ] { y } (12) 

nxl nxm mxl 


where {x} is an approximation of {x} and m < n. 
matrix is taken to be unitary, so that 

The transformation 

[V] T [V] * [I]. 

(13) 

The reduced mth order eigenproblem is then 


[A] {y} * X{y> 

(14) 

where 


[A] - [VJ T [B][V] 

(15) 


and X is an approximation of the eigenvalue X. 

The essence of the reduction scheme lies in the choice of the 
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transformation matrix [V]. In the tridiagonal reduction method, the 
Lanczos algorithm (refs. 11 and 13) is used to build the [V] matrix, 
vector by vector, i.e., 

[V] * [{v^ {v 2 > {v^}] (16) 

such that the reduced mxm matrix [A] , is tridiagonal and its eigenvalues 
approximate the higher end of the eigenspectrum of [B] (or, equivalently, 
the lover natural frequencies of the structure). 

The algorithm yields 


tA] 


a ll d 2 

N N 

N N 

. N N J 

d 2 . a 22 d 3 

X N X 

d 3 a 33 n 

\ N N 

\ X 

n N m 

\ X 

N d x a 

m mu 


X 


(17) 


where the matrix coefficients are theoretically given by the recurrence 
formulas 

» il “ (v 1 } T [B]{v i ) 

; i »l t m (18) 

(v 1+1 ) - [B]{v 1 }-a ii {v 1 }-d 1 {v i _ 1 } 

{v i+l } ■ {; i41 >/tt; i+l ,T{; i+l ,]1/2 

The process is Initialized by choosing an initial trial vector, {v^} 
and setting {v_^} ■ {0}; d^ * 0. 

The Initial trial vector should contain all components of the system 
eigenvectors and must be constrained to eliminate spurious elgensolutions 
(ur •> « ) due to mass matrix singularities. These requirements are 
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satisfied by setting 


{v^ = [B]{w}/[ (Bw) T (Bw) ] 1/2 , (19) 

where (w) is an n-element vector obtained from a random number generator 
routine. 


Reorthogonalization of the Trial Vectors 


Although the trial vectors {v^} generated in equations (18) form 

a theoretically orthogonal set, it has been shown (ref. 16) that they 
rapidly degrade as the computations proceed, such that the later vectors 
are far removed from orthogonality to the earlier ones. This is caused 
by unavoidable computational round-off, which, because of repeated 
multiplications by the unreduced eigenmatrix, [R], tends to amplify 
the contributions of the lower frequency eigenvector components. To 
correct this problem, Gregory (ref. 12) experimented with the use of 
higher precision arithmetic, but found only marginal improvements in 
the final results. Later, Lanczos suggested a single reorthogonalization 
of the trial vectors. While this improves matters substantially, it still 
does not eliminate the precision problem adequately. However, Ojalvo 
and Newman (ref. 14) found that the introduction of an iterative 
reorthogonalization loop can make the trial vectors as orthogonal as 
necessary for extremely large systems. The procedure is as follows: 

Denote each new vector obtained from equations (18) as {v° + ^} and 


iterate, 

{v i+i> ■ <v ° + i’ - j J 1 t(v i ,T(v i+i ,Kv J > 

<vi> ■ ''i+i 5 'A ii v T<v i+i m v 

j*l J J 

until an acceptable vector 

O ■ <vi> - i L t { vX + i>Ky 

j-i J 


is found which satisfies the orthogonality criterion 

max |{v.} T (v® > < lO 2 "* 

l<j<i| J A 

where t Is the total number of decimal digits carried by the computer. 


( 20 ) 
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A normalized form of the reorthogonalized trial vector is finally 
obtained through 


(v m ) ’ 


r s4*l i j r r s+l 1 T / s+1-, 1/2 


( 21 ) 


Experiences gained through application of the FEER program to 
a large variety of problem types and sizes have indicated that an 
average of only two reorthogonalizations are required per trial vector 
generation. 


Size Criteria for the Reduced Eigenproblem 


As a result of numerical experiments and applications (refs. 14, 
15, 17-19), it has been found that in cases where m<<r (where r is the 
total number of structural modes, including rigid body modes, and m is 
the size of the reduced eigenvalue problem), a first grouping of more 
than m/2 lower frequencies of the reduced system are in accurate 
agreement with the corresponding number of exact frequencies, provided 
that m i 7, i.e., when at least seven trial vectors are chosen. The 

remaining reduced system frequencies are spread across the remaining 
exact spectrum, with the last one representing a lower bound on the 
highest exact frequency of the unreduced problem. 

Thus, if the user requests q lower frequencies of the structure, 
the order of the reduced eigenvalue problem is 


r Bin [2q+l,r] ; q > 3 i 
t min [7, r] ; q £ 3 * 


( 22 ) 


It should be noted that in all cases m _< r, and whenever m is 
set equal to r, all the structural modes of the unreduced problem are 
generated. 


Error Bounds on the Computed Eigenvalues 


One of the inherently striking features of the tridiagonal reduction 
method is that the solution of the reduced, tridiagonal eigenproblem 


- X(y> 




[A] {y > 


(23) 



and the off-diagonal elements of [A] automatically provide accurate 
error-bound parameterr for the extracted eigenvalues. In particular, 
it can be shown (ref. 20) that absolute error bounds for each approximate 
root, u >2, are found from the inequality 


-2 . 
u)^ + a 

2 . 

+ a 


- 1 


d mtl 


(24) 


2 

where is an exact system root, d^^ is the (m+l)th off-diagonal 

element of an [A] matrix of order m+1, and y ^ is the last element of the 
eigenvector corresponding to . 


Program FEER Flow Diagrams and Sample Output 


The overall flow diagram for implementation of the tridiagonal 
reduction method in FEER is shown in figure 1. The reduced system 
eigenvalue problem is solved in block 7 by means of a Q-R algorithm which 
takes advantage of the. symmetrical, tridiagonal form of the eigenmatrix 
and the physical modal vectors and frequencies are finally computed in 
block 9. The details of block 6, "Execute Tridiagonal Reduction 
Algorithm", are given in figure2. Block 6.4 and the associated 
peripheral test conditions are used to generate re-start vectors when- 
ever premature vanishing of a trial vector occurs. This is usually due 
to the existence of multiple or near-multiple eigenvalues, as described 
in reference 13. Figure 3 shows a representative eigenvalue table 
produced by FEER. In this example, the order of the stiffness matrix 
was 3,072, while the size of the reduced problem was 41. As shown by 
the error bound listing, FEER generated 21 lower frequencies to within 
an accuracy of .01%, using only 362 seconds of CPU time on a CDC-6600 
computer . 


COMPUTATIONAL EFFICIENCY STUDIES 


A count of the major arithmetic operations expended in FEER is 
summarized in Table 1, where n denotes the size of the stiffness matrix 
in the analysis set, b and b are average semi -band widths of the stiffness 
and mass matrices, respectively, and q is the number of accurate modes 
requested by the user. Each operation is assumed to consist of a 
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multiplication followed by an addition. 

It should be noted that the average bandwidth parameters are used 
primarily to provide a measure of the number of non-zero matrix elements. 
In actuality, FEER employs efficient packing routines which do not 
require a uniform band structure for efficient computational operations. 
It can be seen that the major computational effort involves decomposition 
of the modified stiffness matrix (step 3) and provides the leading term 
of 1/2 nb^ in the total operation count. One of the positive features 
of the tridiagonal reduction method is that only one such decomposition 
is performed regardless of the number of roots required. 

Operation count and storage requirements for several alternate 
eigensolution methods are compared with FEER in Table 2. The purpose 
of this comparison is to provide an indication of the potential 
efficiency of each method, assuming that an equally adept and kncwledg- 
able programmer has had a chance to employ the same time-saving tricks 
in each case. For this reason, several excellent solution techniques 
which achieve high efficiency through special data handling and storage 
methods (see for example, refs. 21 and 22), but nevertheless show a 
high minimum operation count, have not been included in the comparison. 

As in Table 1, the counts are presented in terms of average bandwidths 
which are again to be interpreted as a measure of non-zero matrix entries 
rather than in terms of a specific band structure. 

2 a 

It can be seen that in the Givens method the operation count ( -jn ) 

and the storage requirements as well ( O(n^) ) become prohibitively 
large when the size of analysis set grows beyond more than a few hundred 
degrees-of-freedom. 

2 

The leading term In the Inverse Power Method (NASTRAN) is qnb / 2 
as compared to nb^/2 for FEER, since at least one shift per extracted 
root and a subsequent triangular decomposition is typically required in 
the former method. Based on this assumption and the additional supposi- 
tion that an average of seven Iterations per eigenvector are required 
in the Inverse Power method, theoretical operation-count ratios (Inverse 
Power Method/FEER) are presented as a function of semi-band width and 
the number of required roots in figures 4 and 5 for the cases of 
diagonal and consistent mess matrices. These curves provide only an 
approximate estimate of the relative time savings actually accrued for 
several reasons. First, the structure of the stiffness matrix influences 
the decomposition strategy employed in NASTRAN via the active column 
approach. In addition, there is no a-priori knowledge of the actual 
number of shifts and Iterations which will be required in the Inverse 
Power method for any given problem application. In general, both the 
number of shifts and iterations tend to increase with the number of 
roots extracted, so that the curves indicating improved efficiency 
of the Inverse Power method for a very large number of extracted roots 
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and small bandwldths are unrealistic. 

Table 2 also shows approximate operation counts and storage 
requirements for Gupta's Sturm Sequence method (ref. 23) and a current 
version of the subspace or Block-Stodola method (ref. 9). The storage 
requirements for each of these methods, as well as the Inverse Power 
and PEER methods, are all on the same order of magnitude. In Gupta's 
method the count of 25nb2q is based on his assumption that approximately 
2nb^ operations are involved in examining the Sturm sequence for one 
trial root value, and that about twelve such values must be examined 
for each accurately predicted root (ref. 9). With regard to the 
Subspace Iteration method, the leading term in the count, nb^, is twice 
as large as in FEER and all other terms involving the same functional 
forms of the parameters n, b, E, q are 'also much larger. In addition, 
the reduced eigenproblem which is solved for improvement of the sub- 
space is not tridiagonal so that the count for this operation is on 
the order of q^ as compared to q2 for the tridiagonal reduction method. 
Finally, the assumption of eight subspace iterations may not be very 
reliable, since this depends on the choice of the starting subspace, 
which is somewhat arbitrary. 

Table 3 presents a set of actual computer runs comparing the 
CPU execution times of FEER vs. the Inverse Power and Givens methods 
in the NASTRAN READ module. The results Indicate that the more efficient 
decomposition operations and shift strategy incorporated into Level 15 
have yielded significant improvements in the Inverse Power method as 
compared to the Level 12 version (see also ref. 24). 

However, the run times for comparable or identical problems are 
generally 5 to 20 times faster with FEER than with the Level 15 Inverse 
Power method when between 5 and 20 accurate modes are requested. This 
result is in rough agreement with the operation count ratios shown in 
figures 4 and 5. In problem No. 2, which is relatively small and could 
therefore be treated with the Givens method, the execution time via 
FEER was approximately 3 times as fast, since only 35 modes were re- 
quested, while in the Givens method the user has no choice and must 
pay the penalty of having all the eigenvalues calculated (in this 
particular case, 105). 


CONCLUDING REMARKS 

Significant computational efficiencies are achieved in the FEER 
program primarily due to the tridiagonal reduction method of modal 
extraction. Basically, the subspace of trial vectors generated via 
this method are sufficiently rich in the lower modes to provide a 
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single, reduced, tridiagonal eigenproblem whose solution provides these 
inodes with a high degree or accuracy. This feature distinguishes it 
from the usual subspace or block iteration methods, where the trial 
vector subspace is established somewhat arbitrarily and subsequently 
improved through repeated solutions of reduced eigenproblems. The 
tridiagonal reduction method employs only a single, intitial shift of 
eigenvalues and hence requires only one matrix decomposition. It is 
consequently much more efficient than the Inverse Power Method with 
shifts when more than one or two lower modes are requited. FEER is 
also extremely efficient for out-of-core operations and requires only 
(15,000 + 7.n) central memory words, where n is the order of the 
analysis set. Another feature of the method is that the reduced 
problem is generated automatically, starting with a random trial 
vector, and this avoids one of the basic weaknesses of techniques 
requiring either a judicious selection of starting vectors or retained 
nodes . 


It is concluded that the tridiagonal reduction method used in 
FEER would serve as a valuable addition to NASTRAN for increased 
efficiency in obtaining structural vibration modes. 
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Figure 1* * Overall flow diagram for tridiagonal reduction method 
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Figure 2, - Flow diagram for block 6 . execute tridiaronal reduction 
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